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Dissipation and velocity distribution at the shear-driven jamming transition 
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We investigate energy dissipation and the distribution of particle velocities at the jamming transi¬ 
tion for overdamped shear-driven frictionless disks in two dimensions at zero temperature. We find 
that the dissipation is cansed by the fastest particles and that the fraction of particles responsible 
for the dissipation decreases towards zero as jamming is approached. These particles belong to an 
algebraic tail of the velocity distribution that approaches ~ v~^ as jamming is approached. We 
further Hnd that different measures of the velocity diverge differently, which means that concepts 
like “typical velocity” may no longer be used—a finding that should have implications for analytical 
approaches to shear-driven jamming. 

PACS numbers: 64.60.-i, 64.70.Q-, 45.70.-n 


The hypothesis that the slowing down of the dynamics 
in systems as different as supercooled liquids, granular 
materials, colloids, foams, and emulsions, have a com¬ 
mon origin in the properties of a critical point, point J[l], 
has inspired a great amount of work on jamming the last 
decade. Several models have been used to try and pin¬ 
point the properties of this jamming transition. Some of 
them have centered around a greatly simplified numeri¬ 
cal model of spheres with contact-only interaction. One 
important branch has been to examine the properties of 
randomly generated static packings [2] whereas another 
has been to study the jamming transition through simu¬ 
lations of elastic particles under steady shear [3]. 

A key feature of jamming is the approach of the con¬ 
tact number z to the isostatic number Zjso which is just 
enough for mechanical stability. It has recently been 
shown [4] that this is directly linked to the divergence 
of rjp = p/'j —the pressure equivalent of the shear vis¬ 
cosity. A related phenomenon is the increase in particle 
velocity as (/> —>• 6]. This is related to the distri¬ 

bution of particle displacements due to a small shear 
increment which has been determined both in experi¬ 
ments of sheared granular materials [7] and in quasistatic 
simulations [5, 7]. It was there found that this distribu¬ 
tion is sufficiently wide that the non-Gaussian parameter 
/3 ((Aj/)^) — 1 diverges as (j)j is approached from 
below, granted that the shear step is sufficiently small. 

In this Letter we show that there is more to the par¬ 
ticle velocity distribution than has so far been realized. 
Dissipation is mainly caused by the fastest particles and 
we find that the fraction of particles that are responsible 
for the dissipation decreases towards zero as jamming is 
approached. This behavior is related to an algebraic tail, 
P{v) ~ v~^, in the velocity distribution and we show that 
the velocity histograms determined at the jamming den¬ 
sity approach this limiting behavior as 7 —> 0. Since rhe¬ 
ology and dissipation are linked through power balance, 
the understanding of this phenomenon is right at the cen¬ 
ter of the phenomenon of shear-driven jamming. We also 
note that this finding has profound consequence for ana¬ 


lytical approaches to jamming since it implies that differ¬ 
ent measures of the velocity behave differently and that 
concepts like “typical velocity” therefore become useless. 

Following O’Hern et al.[2] we use a simple model of 
bi-disperse frictionless soft disks in two dimensions with 
equal numbers of disks with two different radii in the ra¬ 
tio 1.4. Length is measured in units of the diameter of 
the small particles, ds- We use Lees-Edwards bound¬ 
ary conditions [8] to introduce a time-dependent shear 
strain 7 = t'j. With periodic boundary conditions on 
the coordinates Xi and j/i in an L x L system, the po¬ 
sition of particle i in a box with strain 7 is defined as 
Tz = {xi + UVijyi)- The ordinary velocity is v-°* = r^, 
but in the following we consider the non-afSne velocity, 
Vi = — VR(ri) where VR(ri) = 'jyiX is a uniform 

shear velocity. With the distance between the centers 
of two particles and dij the sum of their radii, the relative 
overlap is Sij = 1 — Vij/dij and the interaction between 
overlapping particles is V{rij) = e6fjf2-, we take e = 1. 
The force on particle i from particle j is ifj = —ViV (r^j. 
The simulations are performed at zero temperature. 

We consider two different models for the energy dissi¬ 
pation. In both cases the interaction force is ff' = J2j 
where the sum extends over all particles j in contact with 
i, and the equation of motion is 

+ = ( 1 ) 

Most of our simulations have been done with the RDq 
( reservoir dissipation) model with the dissipating force 

fRD.z = -fcdVi. (2) 

We take kd = 1, = 0, and the time unit tq = dgkd/e- 

We simulate N = 65536 particles with shear rates down 
to 7 = 10“^°. Checking for finite size effects at 7 = 10“® 
we found no difference when using instead N = 262144. 
The equations of motion were integrated with the Heuns 
method with time step At = 0.2tq. 

Some additional simulations have also been done with 
the CDo model (CD for “contact dissipation”) with dissi- 
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FIG. 1. (Color online) Velocity distribution and dissipation. 
Panel (a) shows the velocity distribution function Pv{v) vs 
n /7 with simulation parameters {4>,'y) and symbols as shown 
by the inset. The dashed line has slope —3. Panel (b) shows 
the part of the dissipated power which is dissipated by the 
fraction x of the fastest particles. Panel (c) shows the same 
quantity for the CDq model. These data are at 0 = 0.8433 « 
4>j and five different shear rates. The insets of panels (b) and 
(c) are * 50 —the fraction of particles needed to dissipate 50% 
of the power. The inset of panel (b) shows that xzo for the 
RDo model decreases as 4> increases whereas the inset of panel 
(c) shows 2:50 at (pj decreasing with 7 for both RDq and CDq. 

pation due to velocity differences of disks in contact[3, 9]. 
Details of these simulations may be found elsewhere[10]. 

A key quantity in the present letter is the energy dissi¬ 
pation. We here just remark that this is a central quan¬ 
tity due to the relation between dissipation and rheol¬ 
ogy from power balance, Va'p = kd ^.nd we 

therefore believe that the considerations here may be in¬ 
strumental in developing a better understanding of shear- 
driven jamming. 


Our first key result is that most of the energy is dis¬ 
sipated by a small fraction of fast particles and, further¬ 
more, that the fraction of particles needed to dissipate 
a given part of the power decreases as jamming is ap¬ 
proached. Note that “fast” is here used in a relative 
sense. For low 7 all particles are slow, it is only u /7 that 
can be big. To study the dissipation we introduce the 
velocity distribution function Pv{v) such that Pv{v)dv is 
the fraction of particles with velocity u < |v| < u -I- dv. 
Fig. 1(a) shows Pv{v) vs n/7 both at hve densities below 
(j)j, and for three different shear rates at (pj. (To get 
histograms of good quality down to small P^ we use bins 
that are equally spaced in Inn.) The different simulation 
parameters {(p, 7 ) and their corresponding symbols are 
shown in the inset of panel (a). The points connected by 
solid lines and dash lines, respectively, show two differ¬ 
ent ways to approach jamming. The solid line connects 
{(p, 7 ) at ^ < (pj and at sufficiently low 7 to be very close 
to the hard disk limit. The dashed line connects three 
points at 0 = 0.8433 ~ (pj. Here jamming is approached 
as 7 —^ 0 . 

To study the dissipation with focus on the fast particles 
we define 

poc poo 

x{v) = / Pv{v')dv', C2{v) = / Pv{v')v'^ dv'. 

J V j V 

( 3 ) 

Here x{v) is the fraction of fast particles with |v| > v 
and kdC2(v) is the dissipating power due to the same 
particles. We also define C2{x) = C2{v{x)), where v{x) 
is the inverse of x(v). Fig. 1(b) shows the normalized 
C2 vs x for the data in panel (a). The faster particles 
always dominate the dissipation but this effect becomes 
more pronounced—the curves get steeper—as jamming is 
approached; a smaller fraction of particles is then needed 
for a given part of the dissipation. As a simple quanti¬ 
tative measure we introduce X 50 , shown in the inset of 
panel (b), as the fraction of the fastest particles that dis¬ 
sipates 50% of the power. For the hard disk limit (solid 
line) 2:50 decreases as (p increases towards (pj. The be¬ 
havior of X50 at (p = 0.8433 (pj is shown by the open 
squares in the inset of panel (c); 2:50 decreases with de¬ 
creasing 7 and gets as low as 0.16% at the lowest shear 
rate, 7 = 10“^°. We believe that this localization of the 
dissipation to a few faster particles is related to plastic 
events or avalanches that are found above (pj, as already 
speculated by others [5]. 

Panel (c) shows that the CDq model behaves similarly. 
In this model it is the velocity differences of contacting 
particles that is the quantity of interest rather than the 
non-affine velocity, and C2 is defined analogously. The 
main data in panel (c) is C2 (x) at (pj for the CDq model 
which is very similar to the three data sets at cpj in panel 
(a). As a more detailed comparison the inset of panel (c) 
shows X50 against 7 at (pj for both the RDq model and 
the CDq model, and it is clear that this fraction decreases 
with decreasing 7 in both models. The effect studied here 
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is thus not just peculiar to the simpler RDq model. [12] 

The evidence from Fig. 1 strongly suggests that 
C 2 {x)/C 2 {^) approaches a step function as (j) ^ (j)j and 
7 —>■ 0 , and this is the main result from the first part of 
this Letter. For C 2 {x)/C2{1) to approach a step function 
the limiting distribution has to have a tail 

P^{v) - (4) 

since that would make C2(v) diverge. We note that ex¬ 
periments on dense granular flows have led to similar 
conclusions [13]. Before turning to more elaborate analy¬ 
ses we note that the dashed line in Fig. 1(a) with slope 
= —3 gives some support for Eq. (4) as the limiting be¬ 
havior at as 7 —>• 0 . 

For the further analysis it is important to understand 
the origin of the wide distribution. We note that the non- 
afhne velocity in the RDq model is related to the sum of 
all (repulsive) contact forces that act on the particle. The 
non-affine velocity of particle j is Close 

to jamming, the forces on most particles almost cancel 
one another out, and the total force is typically very small 
compared to the average force, <gC /®-, as has also been 
noted by others[14]. There are however some particles for 
which the forces don’t balance one another out, and the 
velocity of these particles can then be much larger than 
the average velocity. The wide distribution is thus due to 
the big difference between the individual forces and the 
typical total force. 

A consequence of this picture is that the maximum ve¬ 
locity is bounded by the typical /fj which means that the 
possibly algebraic distribution is cut off by an exponential 
factor where Vc ~ fij/kd p/{kdds). (This also 

suggests Vc/j ~ r]p.) This behavior is seen in Fig. 1(a) 
as the approximately rectilinear (i.e. algebraic) behaviors 
for intermediate values of Pv{v) turn into more rapid de¬ 
cays at higher velocities. One therefore expects the tails 
in the distributions to be described by P{v) ~ 
and this exponential decay becomes a complicating fac¬ 
tor, when one attempts to determine r from P(v). 

Our second key result is that different measures of the 
velocity behave differently. This is important since it 
means that concepts like “typical velocity”—used in var¬ 
ious theoretical approaches—then become useless. Fig. 2 
shows simulation results for hi and hrms, defined through 

hi = (]v|), and h^^, = (v^) , (5) 

with the notation v = v/ 7 . In panel (a) these quanti¬ 
ties are plotted against 7 and are found to diverge al¬ 
gebraically with different exponents: hrms ~ 7 “/^/^^*^ rv. 
^- 0.34 ^ ^-u^/zv ^ ^- 0 . 28 ^ expressions 

follow by taking h = 7“^/^ in ~ 

7&^)[15] with the scaling dimension uq equal to Uy 
for hi and /3/2 for hrms- The latter follows from r] ~ h^ms 
and 7 ^ 0 ) - ((/)j - <?i)“'^.] 

It is instructive to also examine the same quantities 
with data below (j)j, close to the hard disk limit. The 



FIG. 2. (Color online) Two different measures of the velocity. 
Panel (a) shows hrms and hi for (p = 0.8433 ~ (pj vs 7. Panel 
(b) shows the same quantities for data below cpj plotted vs 
distance to jamming; only the points with cpj — (pes < 0.006 
were used for the fits, (pgfi is the effective density. The lines 
are {<pj — (pes)~^'’ and ~ {(pj — (pe«)~^^'^- Panel (c) connects 
back to Fig. 1(a) but is C 2 , related to hrmsi against Ci which 
is related to hi . From the crossings of the dashed line one may 
read off the relative contribution to hi from the fraction of the 
fastest particles that disspate 50% of the energy—a quantity 
that decreases with decreasing 7. 

starting point is the relations for hard disks, ^[“^(((i) ~ 
{(pj - and h“,((^) -- {(pj - (p)~^/'^, that follow by 

using 6 = 7 in the scaling expressions and considering 
7 —0. In Fig. 2(b) we make use of the effective-density 
mapping of soft disks onto hard disks, 0^'^((^eff) = 
0{(p, 7), where the effective density is (peS = <P — 
with c = 1.53 and y = 1.09, as detailed in Ref. [16]. 
Fig. 2(b) shows hi and hrms against (pj — ^eff- The solid 
line gives the exponent = 1.10 in agreement with 
{'P’J ~ for the particle “velocity” in Ref. [5]. 

The dashed line gives /3/2 = 1.25. (The value fp = 2.50, 
is somewhat low in comparison to recent estimates[I5], 
but this could be due to not including corrections to 
scaling[I5].) Note that the exponents from Fig. 2(a) and 
(b) are consistent when using 1/zv = 0.26[I5]. 

The reason for the different behavior of hi and hrms is 
that the dominant contribution to these quantities come 
from different velocity intervals. This is illustrated in 
Fig. 2(c) which shows how C2{x) in Eq. (3) and Ci{x) 
(for instead of u^) increase to their respective limits 
C 2 (I) = hrms rind ( 71 ( 1 ) = hi, as a; (the fraction of parti¬ 
cles included in the calculations) increases. The different 
curves get steeper for smaller 7 and for 7 = I0“^° we 
find Cl (a;) w 0.05 when C2{x) = 0.5 which thus shows 
that h^j^g gets a considerably larger contribution from the 
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FIG. 3. (Color online) Velocity distribution and the determi¬ 
nation of the exponent r for different 7 at 0 = 0.8433 « <)>j. 
Panel (a) is the velocity distribution P(vy). Panel (b) shows 
Vc and r —3 vs 7 from fitting to Eq. ( 6 ). Fitting r — 3 ~ 7 '^'' for 
10“® < 7 < 5 X 10“^ gives qr = 0.09 ± 0.02. Panel (c) shows 
that rip and Vc are proportional—the dashed line through Vc 
is rip/6.1. Also shown are Vcr —which is where the histogram 
crosses over to exponent —r, and W 50 above which 50% of the 
dissipation takes place. These similarity of these two quan¬ 
tities suggests that the dominant part of the dissipation is 
given by particles governed by the algebraic tail. 


highest velocity part of the histogram than vi. An ex¬ 
trapolation of these curves to the 7 —s- 0 limit would give 
a step function (though this is not as clearly suggested by 
the data as in Fig. 1) which would imply that u^nis and vi 
were controlled by different velocity intervals, and that 
there is no reason for these quantities to be at all related. 

We will now relate our two key results of Figs. 1 and 
2 to properties of the velocity distribution function with 
the goal (1) to examine how the exponent in Eq. (4) ap¬ 
proaches —3 as 7 0 (this exponent will be denoted by 

—r) and ( 2 ) to shed some more light on the mechanism 
that allows vi and Urms to diverge differently. We have 
then found it convenient to use P{vy )—the distribution 
of the absolute value of the y-component. This quan¬ 
tity differs from Pv in that it approaches a constant at 
small velocities—a feature that makes it easier to find an 
analytical expression that fits the data. Fig. 3(a) shows 
P{vy) for several different 7 at ^ = 0.8433 « (/)j together 
with the solid lines that are fits to the expression 


with A, Vc^ Vai Vsi and r as free parameters. This ex¬ 
pression crosses over from a constant at small Vy to a 
large-Uy tail with (as discussed above) and 

the crossover is governed by an additional term in the 
denominator, {vy/vaY- For best possible fits, a should 
be an additional free parameter, but since a anyway 
tends to be close to 2 and a = 2 opens up for analyt¬ 
ical calculations[17], we here fix a = 2 . 

Fig. 3(b) shows the exponent as r — 3 vs 7 . The recti¬ 
linear behavior suggests an algebraic decay, r — 3 7 '?'’ 

with Qr = 0.09 ±0.02, consistent with the limiting behav¬ 
ior of Eq. (4). The points on top of Eig. 3(c) are the cutoff 
velocity, Vc (solid circles) from the fits and rjp (squares) 
directly from simulations. The dashed line through the 
solid circles, which is 77 ^/ 6 .! and not a fit to the data, con¬ 
firms the expectation that these quantities should behave 
the same. The covariation of r and Vc (compare panels 
(b) and (c)) makes the fitting difficult—a small decrease 
in r can be compensated by a small decrease in Vc since a 
smaller r gives a slower decay, while a smaller Vc gives a 
faster decay. This effect is most problematic at the low¬ 
est shear rate, 7 = 10 “^°, and this point is therefore not 
included in the determination of gr¬ 
it is now interesting to determine the size of the re¬ 
gion governed by the algebraic decay, Vy'", and we there¬ 
fore calculate the crossover velocity Ucr = 
that describes the crossover from exponent —2 to expo¬ 
nent —r, by equating the two velocity-dependent terms 
in the denominator of Eq. ( 6 ). Fig. 3(c) shows that Ucr 
behaves about the same as v^o —related to X 50 above— 
which is the velocity above which 50% of the dissipation 
takes place. We thus find that the dissipation is largely 
governed by the particles in the algebraic tail. Recalling 
the conclusions from Fig. 1, it is clear that the fraction 
of particles in the algebraic tail decreases with 7 and 
this is also shown by the big open circles in Fig. 3(a) 
which are P('Dcr) vs Vcr for 7 = 10 “® and 10 “^°—the last 
two points in a persistent trend to smaller P(vct)- The 
fact that the fraction of particles in the algebraic tail de¬ 
creases with decreasing 7 means that their contribution 
to decreases (ui is instead dominated by the slower 
particles) whereas they always give the dominant contri¬ 
bution to u^jjjg[17]. This explains the different behavior 
of vi and firms- We finally note that the algebraic tail 
actually becomes wider as 7 decreases. Fig. 3(c) shows 
that Vc increases faster than Vcr which means that the 
width of the algebraic tail—the region between her and 
Vc — increases with decreasing 7 . 

To summarize, we have found that the fraction of par¬ 
ticles that are responsible for the energy dissipation de¬ 
creases towards zero as jamming is approached. These 
particles belong to a tail in the velocity distribution 
that approaches P{v) v~^ at jamming. We further 
find that different measures of the velocity diverge differ¬ 
ently which means that concepts like “typical velocity” 
no longer appear to be useful—a result of importance for 
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analytical approaches to shear-driven jamming. 
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Supplemental material 

Derivations related to “Dissipation and velocity distribution at 
the shear-driven jamming transition”, Peter Olsson. 


Summary 

The purpose of the calculations below is to illustrate the mechanism that 
gives different behaviors for hi and v^ras- The idea is here to split the expres¬ 
sion for the velocity distribution 


P{Vy) = 


y-Vvlv„ 


1 + {Vy/VaY + [Vy/VsY' 

into three different parts 

1/[1 + iVy/Va)\ 

PiVy) = { l/{vy/Vs)Y 

0 , 


Vqx ^ ^ 


(7) 


Vc < Vy. 


We may then determine hi = {vy) and h^^^g = (vyj by calculating some 
integrals analytically, as shown below. The results are then written in terms 
of the different velocity scales, h^, h^, her, and he- Fig. [Hand Fig. 3(c) show 
that all these quantities to decent approximations diverge algebraically for 
small 7 . We have ha ~ 7 ”'^“, Vs ~ her ~ and he ~ 7 “'^'' with the 

exponents qa ~ 0.18, ~ 0.29, gcr ~ 0.46, and Qc ~ 0.71. 

The hrst moment of the velocity becomes 


Vi = {Vy) ^ Va 


In 


+ 


0.621ogio(l/7) + 


For the second moment of the velocity we get different results when r is 
close to 3. For r not too close to 3 (which is the case relevant for comparisions 
with our simulation data) we hnd 


^rms ~ ^’aVer 

whereas for r —)■ 3 we instead hnd 


1 + 


r — 3 


( 8 ) 


V ~ 

^rms 


'Fa'Fcr “1“ 




VaVer [1 + (Qc “ ?cr) ln(l/7)] 


VaVer [1 + 0.58 logio(l/7)] • 


(9) 


1 













Figure 4: Velocity scales Va, Vg, and Vcr of Eq. ( 6 ). All these quantities behave 
algebraically for the smaller 7 , to decent approximations. 


The small-7 limit 

With the values for the exponents given above, the leading small -7 behaviors 
become ~ VaVcr ~ ^-('?cr+(?a )/2 ^ ^-o.64 very good agreement with 
hrms ~ in Fig. 2(a). 

For hi, the full expression hi ~ ha[ln(hcr/ha) -|- l/(r — 2)] reproduces the 
exponent —0.28 in Fig. 2(a). In the 7 —)■ 0 limit, however, the second term 
vanishes and the slowly changing logarithm may be replaced by a constant, 
and this gives hi ~ 7 “*^“ ~ which is clearly different from the mea¬ 

sured exponent, —0.28. We consider the exponent obtained directly from the 
measured data to be more reliable as it is from an excellent fit to data across 
four orders of magnitude in 7 . In contrast, the “analytical” value for the ex¬ 
ponent (—0.18) is obtained by the dangerous process of extrapolating results 
from fitting to an expression (Eq. ( 6 )) that is used without any theoretical 
justification. 
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Derivations 

To determine vi = {vy) and 'u^ms = need the integrals 


Jo = y P{vy) dVy, h = J Pivy)vy dVy, h = J Pivy)Vy du 
which with Eq. (jT]) becomes 




Ip = 


J P{Vy)vPdVy = 


fVcr 


'0 1 + {Vy/VaY 


rvc ^ ^ ^ 

dVy+ VliVy/VsydVy 

J Vrr 


It is then convenient to consider the two terms above separately. 


First term, 

We here use x = Vy/va, dx = dvy/va, and Xcr = hcr/w and handle the 
different integrals separately for different p: 

• With p = 0: 

(ll r'^cr \ pXcT ][ 

pQ = .A vy = Va ——= ha[arctana:]o"" 

Jo l + {vy/vay Jo 1 + x^ 

= haarctan(hcr/ha) (10) 


• With p = 1: 



I T— 

0 1 + {Vy/Vay ' 


• With p = 2: 




_ y. _ (ji) 

l + iVy/v.y ^ 


y ln(l + {v^jvaf). 


1 1 ^cr 

-ln(l + x2) 

2, Jo 

( 11 ) 


_0 X _0 

/ - -dx = vPx — arctanxL'"'' 

Jo 1 + 

v^Vcr - vl arctan(hcr/ha)- (12) 
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Second term, 

We here get an expression for general p: 



This splits into two cases. Ifr — 1—p3>0 (which is always the case in our 
simulations) and Vc S> her (which allows us to skip the second term) we make 
use of vl~‘^ = ^l/vl to get 


r - l-p 7'52~ p -1 

^a^cr 


^ r — l—p\VcrJ r — l—p 

For r — 1—on the other hand, (possible for p = 2 and very close to 
jamming, r —)■ 3), we get 


~ ln( — ) — ln( — ) = In ( — 

p ^ Vi";/ Vi";/ ^ Vi"; 


Taking both terms together 

The normalization becomes 


Jo = Va arctan 


For the first moment we get 


/. = -In 1 + 


In C- + 


1 hi TT ^ 

H ~ —Va 

r — IVrr 2 


r — 2 


r — 2 


The logarithmic correction is ln(hcr/ha) = ln (7 ”?“) rsj (gcr-ga)ln(l/7)- 

The average velocity becomes (using 0.18 and gcr ~ 0.46) 

vi = (vy) = h/Io ^ Va + 0.621ogio(l/7) + • 

L \ Vn / T A\ L T A_ 
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For the second moment and for r not too close to 3, Eq. fll3p leads to 


h = v„u 


a^cv 


1 + 


r — 3 


— nf arctan 


(17) 


and since arctan x < 7i/2, the second term may be neglected and we get 


Vl) = ^rms ~ VaVcr 


1 + 


r — 3 


(18) 


For r —)■ 3 we instead use Eq. flT^ which gives (with qc = 

0.71). 



~ ~ 1 7 '^C \ ~ ~ / 

VaVcr + — hi — ~ VaVcr I + {Qc 
Va \Vcr/ 

VaVcr [I + 0.58logio(l/7)] • 


gcr)ln(l/7)] 

(19) 
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